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We present experimental and numerical results for displacement response functions in packings 
of rigid frictional disks under gravity. The central disk on the bottom layer is shifted upwards 
by a small amount, and the motions of disks above it define the displacement response. Disk 
motions are measured with the help of a still digital camera. The responses so measured provide 
information on the force-force response, that is, the excess force at the bottom produced by a small 
overload in the bulk. We find that, in experiments, the vertical-force response shows a Gaussian- 
like shape, broadening roughly as the square root of distance, as predicted by diffusive theories 
for stress propagation in granulates. However, the diffusion coefficient obtained from a fit of the 
response width is ten times larger than predicted by such theories. Moreover we notice that our 
data is compatible with a crossover to linear broadening at large scales. In numerical simulations 
on similar systems (but without friction), on the other hand, a double-peaked response is found, 
indicating wave-like propagation of stresses. We discuss the main reasons for the different behaviors 
of experimental and model systems, and compare our findings with previous works. 

PACS numbers: 45.70.-n, 45.70.Cc, 83.80.Fg 



I. INTRODUCTION 



Stress distributions in static eranular materials display 
puzzling characteristics Q, 0, Ga that do not quite fit 
into classical elastic descriptions, and have defied at- 
tempts at analytic modeling for some time already. The 
observation of a pressure dip below conical piles, force 
chains, sudden macroscopic changes in stress patterns 
under slight perturbations, and exponential (instead of 
Gaussian) stress distributions, among other phenomena, 
have triggered intense theoretical and experimental 
work. As a result, a multiplicity of stress propagation 
models have been put forward. The g-model [H El 
assumes diffusive behavior for the vertical stress compo- 
nent considered as a scalar quantity, and gives rise to an 
exponential distribution of stresses. Other scalar models 
in turn predict Gaussian Q or power-law distributed Q 
stresses. By postulating a linear relation between stress 
components 0, a wavelike equation 01 is derived for 
stress propagation, the so called OSL model This 
model reproduces the pressure dip [13, and is consistent 
with stresses in silos [l^] • A memory formalism 0] 
contains as special limits the wavelike and diffusive 
behaviors. Furthermore, a recent description in terms 
of scattering force-chains 15J gives rise to wavelike 
propagation on small scales, crossing over to something 
similar to classical elasticity on larger scales. 
Linear elasticity describes the propagation of stresses 
in terms of differential equations of the elliptic type, 
wavelike propagation corresponds to the hyperbolic case, 
while diffusive behavior is the borderline, or parabolic 
case. These three descriptions give rise to very different 



responses |3|, |^] when a small force is applied on a 
localized region on the upper surface of a packing. 
Linear elasticity predicts a bell-shaped response, having 
a width proportional to depth. A diffusive behavior, on 
the other hand, implies that the width of the response 
scales as the square root of depth. Finally, a wavelike 
propagation would be evidenced by a response that is 
maximum on a diffuse annulus of linearly growing radius 
(the "light-cone") in three dimensions, or by a response 
showing two diverging peaks in two dimensions. 
One might expect such differences to be easily resolved 
by properly designed experiments. However, presently 
available experimental results are not conclusive. 
Some small-scale results support the validity of the q- 
model |lfj|. while recent experiments using photoelastic 
techniques p"?| appear to be in conflict with a diffusive 
picture. The memory formalism has been shown to 
reproduce the stress oscillations observed in later ally 
confined packings [l^. Experiments on sand 0, |20| 
show a single-peaked response function whose width 
scales linearly with depth, as predicted by elasticity. 
However the precise shape of the response does not quite 
match [2(J that of a linear elastic medium. It has been 
noted that even systems in the elliptic regime can have 
two peaks in their response functions |2l| . It appears 
at present difficult, based on available experimental 
results, to clearly validate, or disprove, any of the stress 
propagation models that have been proposed. 
Numerical measurements on disordered packings of 
frictionless disks, that respect the property of iso- 
staticity p3, |2^ | , both on-lattice l24l |25| (with contact 
disorder) and off-lattice [2rl l2ll |28| . show two clearly 
distinguishable peaks in their response functions. How- 
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ever in experiments wavelike (or hyperbolic) response 
functions have only been observed in ordered pack- 
ings , and it has been argued that disorder produces 
a crossover to an elliptic description on large scales. 
Notice that the measurement of response functions in 
granular packings poses a subtle experimental challenge. 
A distinction must be made between the response 
to infinitesimal perturbations G l , derived under the 
assumption that the contact network does not change, 
and G' , the response under small but finite perturba- 
tions, i.e. allowing for rearrangements. On disordered 
isostatic packings, G % is singular [2^, 0) Ell and does 
not have a well defined continuum limit. This is so 
because the propagation of stresses is described by 
random multiplicative processes on these systems. In 
practice this implies that G % takes positive as well as 
negative values, whose modulus grows exponentially 
with distance from the point where the perturbation is 
applied. Thus on isostatic disordered packings, any finite 
perturbation, no matter how small, necessarily induces 
contact rearrangements, because a large and negative 
Green function corresponds to a contact that will open 
upon perturbation. This anomalous sensitivity to per- 
turbation is due to isostaticity, and has been suggested 
as being responsible for the tendency of stiff packings 
to reorganize upon perturbation |22||. Based on very 
general physical grounds, one can expect rearrangements 
to fundamentally modify response functions. Thus 
one should expect experimentally measured response 
functions to strongly depend on the magnitude of the 
perturbation whenever isostaticity is satisfied. Possible 
effects of rearrangements are discussed for example in 
Ref. [2^ (See also recent discussions in Refs. [13, EI)- 
Moreover, isostaticity has only been rigorously 
proven [22L |23j for frictionless packin gs. It is known that 
friction gives rise to indeterminacies [32, [33 , and this is 
not compatible with isostaticity. In real packings friction 
is important, and it is at present not clear whether 
isostaticity applies in some restricted sense, or not at 
all. Recent molecular dynamics results 1 3 ll on frictional 
deformable spheres are not compatible with the pac king 
being isostatic, although previous similar studies |35j 
reached different conclusions. 

Interparticle forces in granular packings have been 
previously measured by means of carbon-paper experi- 
ments P1H, |37l |33| , photoelastic techniques 0, |2^ , 
pressure sensors [l!J [2J|, and high-precision bal- 
ances [3S]]. None of these methods is optimal for the 
determination of response functions, a task that requires 
the measurement of small forces over small regions. 
Several of the methods used up to now require large 
forces to be applied to the packing, while others are 
not able to detect forces on single particles but only 
averages over relatively large regions. The use of an 
electronic balance gives very precise results, but requires 
an extremely time consuming scanning of the bottom of 
the packing. 

In this paper we report our first results from a novel 



experimental technique allowing precise measurements 
of response functions in two-dimensional packings. 
Our procedure consists in measuring the displacement- 
displacement response function p2l 123. 125). that is, the 
displacement produced on a given disk by an upward 
displacement of one of the disks on the lowermost la yer . 
For isostatic systems, it has been shown [23, [3 l25| 
that this quantity is exactly equivalent to the stress- 
stress response function. Several recent numerical 
studies make use of this equivalence to measure stress 
responses [H [H [H El III E|, however this is the 
first time that this equivalence between stress-stress and 
displacement-displacement response functions is used in 
experiments. 




FIG. 1: A disk configuration as seen by the camera. The 
camera resolution is 640 x 480 pixels. The approximate size 
of the viewable field is 29 x 22 cm, implying a spatial resolution 
of roughly 0.45 millimeter per pixel. Only around 150 disks 
out of a total of 400 in the container, fit in the camera view- 
field. The average number of layers is between nine and ten. 

In standard response function measurements, a small 
force / is applied on a point on the top surface of a 
packing (defined as the origin of coordinates) and the 
response function G(cc, y) is defined as the excess stress 
induced at {x,y}, divided by /. In our experiments, 
the response function G(x, y) gives the excess stress 
at {0, 0} (the central particle on the lowermost row of 
the packings) produced by a small force / acting at 
{x,y}. If boundary effects can be neglected, these two 
definitions of G(x, y) should give statistically equivalent 
results. In other words, although for a given sample 
these two ways to measure the response give different 
results, after sample averaging one obtains the same 
function, if translation invariance applies. 
Some advantages of the experimental procedure pre- 
sented in this work are as follows. The measurement of 
displacements can be done with much better precision 
and by simpler means than that of forces. Vertical and 
horizontal displacements provide information respec- 
tively on vertical (G y ) and horizontal (G x ) responses. 
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We obtain information on G x (x,y) and G v (x,y) for all 
{x,y} in one measurement. Moreover, our technique 
does not involve the subtraction of two force patterns, 
a procedure which is prone to error, and more so when 
two large force patterns are subtracted to obtain a small 
response. In our experiments, the displacement of each 
particle is a direct measurement of the green function at 
that point. 

This work is organized in the following manner: Sec- 
tion 111 Al contains a description of the experimental 
setup used for the measurement of response functions. 
In Section III Bl our results are presented and analyzed. 
A comparison with numerical results obtained on fric- 
tionless packings is established in Section IIIII Finally, 
our results are summarized and discussed in Section llVI 



II. EXPERIMENTAL SETUP 

A. Experimental device 

The experimental device consisted of a rectangular 
container made of two parallel 615mm by 320mm plex- 
iglass plates. The back and front plates were respec- 
tively 34 and 10 mm thick, and they were separated by 5 
mm. 400 aluminum disks with a thickness of 4 mm were 
confined between the plates. These disks had diameters 
16,17,18 and 19mm (100 of each). The container was 
placed vertically in order to minimize friction effects be- 
tween disks and walls. All disks were painted black, and 
a white circular label was affixed onto each of them in or- 
der to allow for motion measurement using digital means. 
Labels had a slightly smaller radius than the disks they 
were fixed on. 

The disks' positions and movements were recorded us- 
ing a still digital camera with a resolution of 640x480 
pixels. The camera pointed to a rectangle of 29 x 22cm 
around the middle of the plexiglass container, implying a 
resolution of roughly 0.45mm per pixel. Before each ex- 
periment i, the packing was shaken and allowed to settle 
under gravity alone. The disk configuration before per- 
turbation was then recorded using the camera, obtaining 
a pre-image Fig. ^shows a typical image. Next the 

central disk in the lowest row was displaced upwards by 
1mm. This disk was fixed to a micrometric screw that fit- 
ted through a specially devised hole in the bottom border 
of the container. At this point the final configuration of 
disks was recorded and stored to post-image This 
procedure was repeated a total of 330 times. 
A C program processed these images to obtain individ- 
ual disk motions. Our program took pairs of images (be- 
fore and after the perturbation) in b&w (lbpp) format 
as input. For each pre-image, all clusters of white pix- 
els were first identified by the method of burning. Their 
geometric centers were taken to be disks centers. These 
centers were then taken as starting points for the identi- 
fication (burning) of clusters on the corresponding post- 
image. This allowed to find the displacement suffered by 



each disk in the pre-image. For each disk, horizontal and 
vertical displacements G x and G y were calculated from 
the differences between geometric centers. These num- 
bers (two per disk) give the displacement-displacement 
response function for a given packing, which for friction- 
less systems equals the force-force response functions, as 
discussed somewhere else [52, 0> H5- Because of the 
existence of friction in our experiments, displacement re- 
sponses are not exactly equal to stress responses, however 
these differences will be regarded as small and ignored in 
the following. 



a) 




b) 

FIG. 2: Response functions experimentally obtained after 
averaging over 330 packings of 400 disks. The displacement 
displacement response measures the displacement ((a) verti- 
cal, (b) horizontal) of a disk at {x, y} in the bulk, produced 
by the upwards motion of a disk at {0, 0} in the lowermost 
layer. On isostatic systems this is exactly equal to the the 
the vertical excess force on the lowermost central grain at 
{0, 0} when a (respectively vertical or horizontal) unit force 
is applied at {x,y}. 



Response functions were averaged over 330 repetitions of 
the experiment. In order to take averages we subdivide 
the image into square cells, adding values of G only to the 
cell to which the center of the corresponding disk belongs 
before the displacement. 
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B. Experimental results 

Fig. [2^,b show our main results, respectively vertical 
and horizontal response functions. The horizontal re- 
sponse G x equals minus the excess compressive force at 
the bottom produced by a unit force acting in the positive 
x direction in the bulk. Fig. [2b indicates that, on average, 
compressive forces on the bottom increase on the right 
side of the load's application point, and decrease to the 
left of it. The q-model has no prediction for this quantity, 
as it only handles the vertical component of the force. It 
would be interesting to compare these results with the 
predictions of other competing stress-transmission theo- 
ries. 

As Fig. shows, the vertical response G y has the form 
of a bell-shaped curve with a single peak, whose width 
grows with the distance from the perturbation point. 
Two models that predict a single-peaked response are the 
g-model [H 0, 0] , and classical elasticity. Within linear 
elasticity, the width u>(y) of the response grows linearly 
with distance, while in diffusive models like the q-model 
it grows with the square root of distance. 
The width at half- height ui(y) of the vertical response 
function G y can be calculated from the data in Fig. 01. 
Our results are displayed in Fig. Ob- Taking a, 6, c as 
free parameters we fit uj(y) = a * y b + c and obtain 
b = 0.51 ± 0.08. This result, taken at face value, sup- 
ports diffuse behavior of stresses. Assuming diffusive be- 
havior (b = 1/2) and fitting w(y) = {2Dy) 1 / 2 + c, we 
obtain D = 95 ± 5mm. This last result is not in very 
good agreement with the diffusive g-model theory, which 
predicts D « grain size/ 2. Given that we have equal 
numbers of disks of size 16, 17, 18 and 19 mm, the average 
size is « 17.5 mm. Thus we find D « 5x grain size, a 
factor of 10 of from the theoretical prediction. With this 
evidence, we believe that a parabolic per se fit cannot 
be taken as strong evidence in favor of diffusive behavior 
at large scales, given the small number of layers (around 
10) that we have in this experiment |4^. Notice that 
an asymptotically linear behavior of u)(y) is also consis- 
tent with our data (dotted line in OH, if deviations in the 
first few layers are ignored. Our preliminary conclusion is 
then that parabolic widening holds on very short scales, 
however with a possible crossover to linear widening on 
larger scales. 

In Ref. ^(| it is argued that an isostatic system of fric- 
tionless disks should behave according to the predictions 
of the g-model, i.e. diffusively. This expectation is not 
verified in simulations [22I l24l l2fij in which the polydis- 
persity is small and all disk centers are located on the 
sites of a triangular lattice (notice that, although disk 
centers are on a regular lattice, these systems have strong 
contact disorder). For these (on-lattice) isostatic sys- 
tems, the average response function of frictionless hard 
disks shows two peaks that diverge linearly. This behav- 
ior is consistent with theoretical arguments j2^j suggest- 
ing that the assumptions leading to wave-like propaga- 
tion of stresses are exact on isostatic packings. 
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FIG. 3: a) Experimentally measured vertical response func- 
tion G y at different heights y above the bottom. The dashed 
lines are Gaussian fits, b) Width lo of G y (x,y) vs height y 
above the perturbation point (squares). The solid line is a 
parabolic fit w(y) = (2Dy) 1 ^ 2 + c resulting in D = 95 ± 5 
mm. The dotted straight line has slope 0.62, and fits the 
last five points. The inset shows the rescaled response func- 
tion G = LuG(x/u>,y), for all values of y. The solid line is a 
Gaussian fit. 



However it has been argued that the response might be- 
come single-peaked when the disorder is large |29|. In 
order to explore the effect of disorder in the positions of 
the disks, and for the sake of comparison with our own ex- 
periments, we performed numerical simulations to mea- 
sure the response on systems of frictionless disks with the 
same distribution of radii as in the experiments. In these 
simulations, disks do not occupy the sites of a regular 
lattice. This is the subject of next section. 



III. NUMERICAL EXPERIMENTS 

The numerical experiments start by pouring disks, one 
by one, into a rectangular die, following a steepest de- 
scent algorithm 0i EH- The equilibrium position for 
each grain is attained when its center of mass falls be- 
tween the centers of two already deposited grains. This 
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way of packing originates a sequentially deposited iso- 
static structure, however not necessarily a stable one 
(positive stresses)as is the case for the on-lattice sim- 
ulations of Refs . I22L l24l |25| , or the adaptive simulations 
in Refs. |27], l28l . The geometric parameters (size and 
number of disks, size of the container, etc.) used in the 
simulations reproduce the scale of the real experiment 
previously discussed. There is no friction in our simula- 
tions. 

Once the assembly is ready, the central particle at the 
bottom of the die is displaced upwards. Upon perturb- 
ing the system we do not allow for rearrangements, that 
is, we keep the list of contacts unchanged. The isostatic- 
ity of the contact network then allows one to calculate 
the displacement of all other particles ver y s traightfor- 
wardly by upwards propagation |22l l24l l25l |4l| . Because 
of the conservation of contacts our results are relevant for 
the limit of very small perturbations in frictionless sys- 
tems. In experiments, rearrangements are very difficult 
to avoid when perturbing the system as we do, unless 
displacements are exceedingly small. 
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FIG. 4: Numerical results for the vertical (a)) and hori- 
zontal (b)) response function, on model frictionless packings 
with the same geometry as in experiments. These responses 
are calculated assuming the limit of infinitesimally small per- 
turbation, in which contacts between disks do not change. 

Because of the existence of multiplicative effects (2f| in 
the simulated experiment, we find the expected large fluc- 



tuations in the measured response functions, thus aver- 
ages are taken over 10 7 samples. Fig.0]shows the results 
for the response functions, as measured in our simula- 
tions. It is clear that, also in this case as for regular 
packings H2, |2J, |25| G y presents a double-peaked shape, 
characteristic of a hyperbolic 0, |2^, ^3] behavior. Thus, 
the amount of disorder considered in this work does not 
produce a single-peaked response in isostatic systems, al- 
though the equivalent but frictional experimental system 
(Fig. [21 displays a single peak. Comparable results are 
found with a more elaborate but time-consuming adap- 
tive algorithm [27L El l3l] . 



IV. DISCUSSION 

We have measured displacement response functions, 
both experimentally on arrays of frictional disks, and 
numerically for disks without friction. Because of the 
virtual work principle, for isostatic systems the displace- 
ment response equals the stress response. Our simula- 
tions do have this exact symmetry, however in experi- 
ments the existence of friction makes the system not iso- 
static in general j^. Thus the displacement response is 
not necessarily equal to the stress green function, how- 
ever we expect the differences between them to be small. 
In this work we have then assumed that stress responses 
take similar values to what we find experimentally for 
displacement responses. The experimentally measured 
response (Figs. and (SJ has a Gaussian-like shape, and 
its width scales approximately as the square root of the 
depth. This appears as consistent with predictions of 
the g-model However, this model predicts a 

diffusion coefficient D whose value is of order half the 
average particle size. The diffusion coefficient D that we 
obtain from fitting our experimental data (Fig. |3J) is ten 
times larger than this prediction. We notice that, up 
to now, diffusive behavior has been clearly seen only on 
experimental systems of no more than ten layers. We 
must thus remark that our results are also consistent 
with a crossover to linear broadening at large scales. The 
shape of the response is better approximated by a Gaus- 
sian than by a Lorentzian, as isotropic elasticity would 
imply. Recent experiments on sand [Tgl |20| find linear 
broadening of the response on scales of the order of 100 
layers. Similarly in those experiments the precise shape 
of the response is not Lorentzian. Photoelastic exper- 
iments 0, |2^| on somewhat smaller piles also suggest 
linear broadening. 

Clearly, larger systems must be studied before stronger 
conclusions can be drawn from experiments like the one 
reported in this work. Experiments on more extended 
systems are under course at present. However in view 
of the present results, as well as those of previous inves- 
tigations 

EE IH 13 HI HI, a possible scenario is to 
have diffusive behavior at short scales, crossing over to 
some sort of effective linear elastic behavior (with linear 
broadening of response) at larger scales. This would be 
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in line with the expectation that linear elasticity should 
be essentially correct at large enough scales. 
Our numerical results, on the other hand, show wave- like 
propagation of stresses, evidenced by two diverging peaks 
in the response function. There is no sign of crossover to 
a single-peaked response. This is consistent with previous 
results [271 1281 131| . and confirms that, on frictionless iso- 
static systems, disorder does not produce a single-peaked 
response. Our numerical results also show that the sim- 
plifying assumptions leading to the g-model are not jus- 
tifiable for frictionless polydisperse systems. 
We remark that our numerical results are valid for the 
infinitesimal response function, i.e. when contact rear- 
rangements can be ignored (2^, EH EI • In practice the 
limit of infinitcsimally small perturbation may be very 
difficult to attain in experiments, so it would be desir- 
able to have a means to quantify the effect of rearrange- 
ments in numerical simulations of responses. Recent ex- 
periments study the effect of rearrangements in packings 
subject to relatively large perturbation |4^ |. 
On the other hand, in most experimentally feasible se- 



tups, the existence of friction partially removes the iso- 
staticity properties of granular packings |3J| . In the pres- 
ence of friction the system of equations in terms of inter- 
particle forces becomes indeterminate |32t l3^ | and ac- 
cepts a multiplicity of solutions. Thus friction may be 
seen as an additional source of randomness (apart from 
geometric and contact disorder). Numerical considera- 
tion of friction effects normally makes simulations very 
time-consuming (33,IH,0|. Clearly it would be interest- 
ing to include the effect of contact rearrangements, and 
friction, in a realistic but efficient way in simulations. 
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